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Abstract 

We explore the spectrum of a flavor singlet pseudoscalar meson, rf, in two-flavor 
(Nf = 2) lattice Quantum Chromo Dynamics (QCD). The continuum-like relation 
between the topology of the QCD vacuum and the U(1)a anomaly, which prevents the 
rj' meson from being a Nambu-Goldstone boson, is expected to hold in the domain wall 
fermions (DWF) used as a lattice quark field in this work. Although our simulation is 
limited to relatively heavy quark masses and the statistical error is not small despite the 
improvements in the measurements and fitting procedures for meson propagators, we 
obtained m T] i = 819(127) MeV for the Nf = 2 QCD, where the error is only statistical. 
Several sources of systematic errors, which may be significant, are discussed. Results 
for the other mesons are also reported. 
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§1. Introduction 



One of the most fascinating puzzles in the meson mass spectrum is U(1)a problem: 
why the mass of the flavor singlet pseudoscalar meson, rj', is so heavier, m v i = 957.78(14) 
MeV, than that of its flavor nonsinglet counterparts in nature, m n o = 134.9766(6) MeV, 
m K o = 497.648(22) MeV, and m v = 547.51(18) MeV. 1 ) 

Nonsinglet mesons behave as Nambu-Goldstone (NG) bosons with the spontaneous break- 
ing of SU(3)a symmetry in the quark massless limit (m quar k — > 0), ignoring the small QED 
effects whereas 7/ is not an NG boson because U(1)a symmetry is broken by the quantum 
effect, the U(1)a anomaly. The nonvanishing divergence of the flavor singlet axial current, 
A®(x), in the axial Ward-Takahashi identity (AWTI) occurs for an operator O in the case 
of degenerate quarks up to the contact term, 



d^A^Oj = 2m^ lk {P a (x)0) + 5 b>0 2N f { Ptop (x)O), (1-1) 

and expresses the anomalous breaking of chiral symmetry in the last term, which is propor- 
tional to the topological charge density, p t0 p(x). For a sufficiently smooth gauge field, 



Ptop(z) = T^tuvpoteF^Fp^x). (1-2) 



The difference in the pseudoscalar meson masses between the flavor singlet sector, m v >, 
and the non-singlet sector, m n , was estimated by the Witten-Veneziano (WV) relation,, 2 )' 3 ) 

ml = m\, -m\ = ^fxtop (1-3) 

at the limit of N c — > 00. Here Xtop is the susceptibility of the topological charge (Qtop) : 

Ql 



1 top 

Xtop = yj, , Qtop = / Ptop(x) d x, (1-4) 

in pure Yang-Mills (YM) theory in a four dimensional volume VT. A recent result in N c = 3 
YM theory with the overlap fermion 4 ) shows that x = (191(5)MeV) 4 . r]' mass from this 
estimation for Nf = 3 and in the chiral limit, m\ — > 0, is m v t m 970 MeV, which is very 
close to experimental values. 

The direct numerical calculation of the 7/ spectrum is important for checking theoretical 
scenarios such as the WV relation, and should result in its correction in finite N c and nonzero 
quark masses. 

Simulations of rj' physics in pure YM theory with quenched Wilson fermions were carried 
out in pioneer works. 5 )' 6 '' The relation between the topological charge and the mass of 7/ 
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was also explored. 7 ) Unquenched simulations 8 -* 12 ) were performed for two-flavor and for 
2+1-flavor 13 ) of Wilson fermions. Using staggered fermions, has been calculated for 
N f = 0, 2 14 ) and N f = 2 + 1. 47 )> 48 ) Recently there are other interesting investigations, such 
as using twisted- mass quarks 49 ) or a local imaginary ^-term. 50 ) 

In this paper, we discuss the mass of 77' in Nf = 2 QCD with domain wall fermions (DWF). 
is one of the lattice chiral fermions, which has both flavor and chiral symmetries 
even at finite lattice spacing (a > 0), and is thus suitable for investigation of nonperturbative 
physics of chiral anomalies. These features of DWF make their use preferable to the other 
alternative methods of discretization. Wilson fermions break chiral symmetry at a > and 
discretization errors start at 0(aAqcn)- The singlet flavor meson in staggered fermions is a 
very important subject as it may be related to the potential issue about the locality of the 
formalisms in the continuum limit. 16 ' 1 ' 1? ) 

Chiral and flavor symmetry are particularly important for 7/ physics, and the DWF is 
the natural choice of lattice quark in investigations. Chiral symmetry in a DWF is not 
realized perfectly, it is broken due to its finite extent in the fifth direction, L s . The amount 
of breaking can be measured by a shift in quark mass: m quar k = rrif + m res , so that the 
nonsinglet axial current is conserved at m quar k = 0. m res is called the residual quark mass 
and vanishes at large L s for a sufficiently smooth gauge configuration. 21 ) 

Although it is desirable to take L s — > 00 limit, to reduce the computational cost, we re- 
strict ourselves to finite L s = 12 with the combination of DBW2 improved gauge action, 22 )> 23 ) 
which smoothen gauge field at short distance and reduces m res significantly. 24 ) 

The RBC collaboration examined the first large scale dynamical DWF simulation. 25 ) 
Pseudoscalar meson masses and decay constants were computed and fit to the chiral per- 
turbation theory (ChPT) formula. m n , mjc, Tn p , f n , and fx calculated in their work are 
reasonably consistent with values obtained in experiments. The J parameter is closer to the 
phenomenological value than the value obtained in the quenched simulation. The nonsinglet 
scalar meson, do, mass and the decay constant have also been examined both in dynamical 
QCD and partially quenched QCD using partially quenched ChPT. 26 ) 

We will mainly focus on the 77' meson in this paper, but we will also report on the results 
of other mesons belong to other Lorentz and flavor representations, and also investigate the 
signal of the mesons in their excited state and their decay constants. 

As the results are limited to the isospin symmetric case and the number of dynamical 
quarks is two, our focal interest in this paper is to provide a benchmark calculation for the 
study of the general meson spectrum on a dynamical DWF ensemble with various (smeared) 
meson field using a larger statistical sample than in the previous study. 

In §2, the theoretical expectations of 77' meson physics are summarized. We explain the 
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details of the simulation including improvements in the signal-to-noise ratio and the fitting 
methods used to relate the simulation data to physical quantities in §3. The numerical 
results are presented in §4 with a list of their systematic uncertainties. We will summarize 
in §5. 

§2. Theoretical results on physics of flavor singlet meson 

In (continuum Euclidean) QCD with Nf degenerated quarks, the operator of the flavor 
singlet pseudoscalar meson, 77': I(J P ) = 0(0~), is defined by quark operators, qj, as 



where / = 1, Nf is the flavor index. The rf propagator consists of two parts: 



d 3 x(r,'(x,t)r)'\0,0)) = C 75 (t) - N f D l5 (t), (2-2) 

r I 1 Nf ■ * > 

C 75 (t) = -Jd 3 x^—J2 q f (x, t) 75 g / (f,t)g / (0,0) 7 5 g / (0, 0) 

/ Nf ^ Nf ^ ^ 

D, 5 (t)= I d 3 x/^-J2 'ff& *)75?/(^ ^~ X; 0)75? g (0, 0) 
\ * / * a I 

The braces represent the contraction of the quark propagators, S q (0,t). Thus, for example, 
C l5 (t) is (S q (t, 0)75^(0, ^75), the same as the nonsinglet meson (pion) propagator, and 
Dry 5 (t) is the correlation function between disconnected quark loops, which exists in the 
flavor singlet mesons. When D JS (t) is suppressed by the OZI rule, it propagates and acquires 
U(1)a anomaly. 

In dynamical QCD, in which the mass of quark polarizing the gluon, m sea , is equal to that 
of the valence quark consisting the meson operator, m va i> the r/' propagator is an exponential 
function of time with its damping factor being the mass of the meson, m' 

J d 3 x( V \x,tW\6,0)) = C 15 {t) - N f D 15 {t) = A/e-"V + . . . , (2-3) 

at large t. 

A model of the 77' propagator is depicted in Fig. 1. The meson propagator is expressed 
as a series expansion in number of the quark loops with signs reflecting the Grassmannian 
feature of the quark, and the blobs at the ends are the meson operators (2-1). The wavy 
lines connecting the quark loops represent the coupling between disconnected loops attached 
to the pseudoscalar density, which is related to the U(1)a anomaly. 
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Fig. 1. Diagram of rj propagator. 

The meson propagator in momentum space can be calculated from the model. The 
first term in Fig. 1 is the same as the nonsinglet pseudoscalar meson (pion), l/(p 2 + m 2 ), 
and the second term is given by two pion propagators coupled to each other by gluons, 
l/(p 2 + m 2 ) x ml/Nf x l/(p 2 + m^), whose coupling we parameterize as m^/Nj. There are 
Nf combinations of quark loops in the second term. Repeating the similar identification of 
connected pion propagators to the nth order, the momentum space representation of the rj 
propagator can be written as a geometrical series: 

W(p)vH-p)) 

1 1 rrin 1 „ 2 1 mn 1 rrtn 1 



OC -/V f „ Oat O O 



— i v f — r~ l v f ■ — 

p2 _|_ m 2 p2 _|_ m 2 jy^, ^2 _|_ m 2 / ^2 _|_ m 2 jy ^2 _|_ m 2 jy^, ^2 _|_ m 2 



= ^_V _^L 2T> (2-4) 

The m„- pole in the connected diagram, C 75 , is exactly canceled by part of the disconnected 



mi + m 



0) 



diagram, D 75 , and thus the square of the rj meson mass, m 2 ,, is identified by 
which means rj does not behave as an NG boson. In terms of this model, rj spectroscopy 
calculated in lattice simulation should reveal the magnitude of wig, and if it is consistent 
with the WV relation (1-3). 

In a later section, we will also calculate the ratio between D 7B (t) and C 7s (t) ~ A n e~ m7rt . 
From (2-3), the ratio at large t should behave as 

NtD (t) e"V* 4- e ""V( T -*) 

+••• (2-5) 



C 75 (t) e m^t + e -m w {T-t) 

(T ^T 1 - Be-*" 



1 _ Be - Amt + • • ■ , (2-6) 



Am = m v / — m n , B = (2-7) 

The ratio at large t exponentially approaches unity with the exponent being the mass 
difference between rj and pion, which is a signature of the dynamical sea quark. This is in 
contrast to the quenched QCD, in which m sea is taken to be infinitely heavy while m va i is 
kept finite. In this nonunitary theory, the third and higher terms in the quark loop expansion 
(2-3) are missing due to the decoupling of the sea quark, and the resulting meson propagator 
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has an unphysical double pole, 

1 1 m 2 1 

{T}\v)r] {-V)) quenched °C -j- j - ^/^- 2 l^T^ 2 " ( 2 ' 8 ) 

The ratio D^ 5 (t)/C 15 (t) in this case behaves as a linear function of time, 

N D (quenched) , , , 

which is clearly different from (2-6). Thus, to obtain a physical 7/ we must simulate the 
dynamical theory (m sea = m va i). We will examine m v > in the domain wall QCD only at the 
dynamical points. 



53. Simulation details 



3.1. Domain wall fermion (DWF) 
The DWF action is defined as 

S F = ^ ^(x,s)D DWF (x,s;y,s')i/j(y,s'), (3-1) 

x,y,s,s' 

£>dwfO, s; y, s') = 5 s , s '£)I, a + &x, v D^, (3-2) 
X 4 

D L = 2 E K 1 " 7m)^(«) + (1 + -r*KW»-te] + ( M s - 4)5,,,, (3-3) 

£>j> = 2 t( X ~ + (1 + 75)^-1^] 

+ ^ [(1 - 75)<Wi V + (1 + 75)^,0^.-1^] , (3-4) 

where ip(x, s) is a DWF that is located in five dimensional space, (x, s), L s is the size of the 
fifth direction, and the parameter M 5 is the domain wall height. By setting M 5 in a region 
around [0,2], from (3-4), left-(right-)handed zero modes are localized around s = 0(L S — 1) 
and the zero modes undergo exponential damping as s, (L s — 1 — s) increases. When a 
four-dimensional fermion and antifermion, q(x) and q(x), are defined as 

q(x) = 0) + ^p-tj>(x, L s - 1), (3-5) 

q(x) = $(x, 0)^^ + L s - 1)^^, (3-6) 

chiral symmetry is fulfilled even with finite lattice spacing (a > 0) at the L s — * 00 limit. 

However, in the simulation, L s is restricted to be finite, and the AWTI is modified from 
its expression in the continuum theory to, 

djA b ^x)0) = 2{m f + m TCS )(p\x)0) 1 (3-7) 
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i.e., the physical quark mass is shifted to m quar k = rrif + m Ics . m rcs is a small lattice artifact 
called the residual quark mass, defined as 



m ms = hm -, r , (3-8) 



where J§ q (x,t) is an operator similar to the pseudoscalar operator but made of fermions at 
the midpoint of the fifth direction, 20 ) s ~ L s /2, thus the numerator of (3-8) includes the 
contractions between the surface fermions at s = or s = L s and the midpoint fermions 
at s ~ L s /2. For the flavor non-singlet case, b ^ 0, m rcs is an exponential function of L s 
as a consequence of the exponentially localized zero modes to the surface, and vanishes as 
L s — > oo. 

One could further argue 52 ) that the effective Lagrangian contains the diverging, 0(l/a), 
discretization error, which can be corrected by the small shift of the quark mass, m quar k = 
mj + m rcs . The remaining error is 0(a), similar to that of Wilson fermions, however, it is 
an exponentially small number, e aLs , or 0(m rcs ). Although m rcs is small compared with 
the statistical errors we will have in most of observables, we will treat the shifted quark 
mass m quar k = rrif + m res as the physical quark mass so that our analysis is precise modulo 
0(m rcs a, a2), which is a few percent in our simulation. 

On the other hand, for flavor singlet {b = 0) case, J§ q (x,t) in (3-8) can be attached 
to a quark loop that does not propagate in the entire L s in the fifth direction, and is free 
from suppression, the counterparts of m res in the flavor singlet case remains finite even as 
L s — > oo, and reproduces the following anomalous term 51 ) 

J2 (JIM t) ) - <W E (pto P (x, t)0) . (3-9) 

X X 

In summary, DWF even for finite L s correctly reproduces the quantum anomaly of axial 
symmetry with small error due to lattice discretization. 

3.2. Ensemble: actions and parameters 

We employ the Nf = 2 QCD ensemble 25 ^ with DWF actions described in the previous 
subsection. Our gauge action contains an improvement in the sense of the renormalization 
group invar iance, DBW2: 23 ' 1 



Sg= 3 



(3-10) 



i - 8ci) ft^t 1 - P ^( x )] + c iE ReTr ^ ~ R ^(%)] 

P lw (x) = U lt (x)U v (x + ii)Ul(x + i>)Ul(x), (3-11) 
R„u(x) = U^xp^x + p)U v {x + 2p)UUx + fi + u)UUx + u)Ut(x), (3-12) 
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Table I. Lattice ensembles and simulation parameters. 



p 


Cl 


V x T a- 1 [GeV] a [fm] 


Va? [fm 3 ] 


?7l res 


0.80 


-1.4069 


16 3 x 32 1.537(26) 0.1284(22) 


(2.054) 3 


0.00137(4) 


m f 




begin-end(step) traj. 


#config. 


N ■ 

1 ' noise 


0.02 


0.5121(36) 


656-5351(5) 


940 


1 


0.03 


0.5984(31) 


615-6205(10) 


560 


3 


0.04 


0.6415(33) 


625-1765(10), 2075-5615(10) a 


470 


2 



( a ) p or the to/ = 0.04 ensemble, we do not use trajectories 1775-2065 due to a hardware error on trajectory 
1772 that was not detected until lattice generation was finished. 



with (3 = 0.80 and c\ = —1.4069. The parameters of the DWF action (3-4) are set as L s = 12 
and M 5 = 1.8. We measure observables on a 470-940 lattice configuration samples for three 
different masses, rrif =0.02, 0.03, and 0.04, which correspond to m^/rrip pa 0.51-0.64. The 
lattice size is 16 3 x 32, the lattice scale is a -1 ~ 1.5 GeV (a pa 0.13 fm), and the residual 
chiral breaking m res = 0.00137(4) which is about an order of magnitude smaller than the 
input quark masses. Throughout this paper we estimate the statistical error using the blocked 
jackknife method. The size of the block is determined to be 50 trajectories by monitoring the 
autocorrelation of the hadron propagators. A summary of lattice ensembles and parameters 
is given in Table I. Other results on these ensembles can be found in 25), 26), 30)-32). 

3.3. Improvements: smearing and sources 

Before constructing the meson propagators, we describe an improvement for the quark 
propagators in this section. It is known to be difficult to reduce the statistical error of the 
flavor singlet meson spectrum. As we have seen in the previous section, the meson propagator 
includes the correlation function between disconnected loops, D r (r = 75, 7t, 1, 7574, 7i7j): 
whose statistical fluctuation is very large, particularly for large t as we will see. We have 
implemented smearing for a quark operator in a gauge- covariant manner called Wuppertal 
smearing. 27 ^ The smeared quark operator qs is a gauge- covariant superposition of the local 
quark operator q^\ 



Ql&t) - g§(2,t) = ^F c < c '(f,f)g2'(y,t), 



y,c 

2 3 



N' 



i=l 



F c ' c (x,y) = 

[l]x,c;i/,c' & ^S,yi 



[Vi]a?, c; £c' = Ui(x, t) c ' c 8 g+ i s - 5 C ' C Sgtf, 



x,c;y,C 



(3-13) 

(3-14) 

(3-15) 
(3-16) 
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[Vjkc;^ = Uj(y, ty% ,,, - 5^'5^. (3-17) 

The shape of qs in terms of is Gaussian with width u as N — > oo. We set uj = 4.35 and 
N = 40. The overlap between the ground state and the meson operator made of smeared 
quarks is expected to be larger the meson made of unsmeared quarks, and the excited state 
contamination is suppressed for small t, where the statistical error is smaller. 

Both the quark correlation functions, Cp(t) and D r (t), are calculated for a complex Z 2 
noise source, £, defined by 

e n \x,t) = ±[&\2,t)+itPfrt)}, (3-18) 

where n = 1,2, . . . , N noisc are random noise ensembles and £i and £ 2 take values of ±1 
randomly. t) is statistically independent of space-time: thus, it satisfies 

-i noise 

w lim w~ E e (n) (^oe (n) (y,o = o, (3-i9) 

-"noise - > °° iV no i se 

n=l 

-i ^noise 

w lim ixr— E e (n) (^,oe (n) *(y,o = feA*'. ( 3 - 2 °) 

"noiso^ 00 i 'noise -, 
n=l 

which is useful for calculating the disconnected loops as we will see in the next subsection. 
We use the source restricted to a time slice (wall source) for Cr{t) and a space-time volume 
source for Dr{t), and iV noise = 1, 3, and 2 for m/=0.02, 0.03, and 0.04, respectively. 

3.4. Meson operators and correlation functions 

Our naming convention for meson fields is similar to that used by the particle data 
group, 1 - 1 but our simulation is limited to having only up and down quarks (TV/ = 2) with 
degenerate masses and zero electric charges; thus, the meson spectra are inevitably different 
from those in the real world. The systematic error from these omission may be comparable 
or smaller to our target precision of ~ 10 %. This point certainly needs further investigation. 

The Hermitian interpolation fields for flavor nonsinglet meson in our simulation, ir, p, ao, 
a\, and bi, and singlet fields, r]', to, f , f\, and hi are defined in terms of quark operators, 
qjj and qjj as follows: 

1 2 

tt/(^0 = -t= Yl ?/,/(^*) 7 /^ i 769/,ff(z,*), (3-21) 

/,9=1 

1 3 2 

p 7 (f, i) = -=J]^ qij(x, t)r b f Jjiq I>g (x, t) (3-22) 
V i=i f,g=i 
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Table II. Meson operators in the simulation and their quantum numbers. 



Meson type 


jPC 


r 


nonsinglet 


singlet 


pseudoscalar 


o-+ 




7T 


»/ 


vector 


\— 


ili a 


P 




scalar 


0++ 


1 


Oo 


/o 


pseudovector 


1++ 


ilbli a 


ai 


./l 


pseudovector 


1+- 


i-lilj a 


&i 


hi 



( a ) average over i, j = 1, 2, 3 is taken. 



1 2 

a 0I (x,t) = —= ^ qi,f(x,t)T b f>g q I>g (x,t), 



ai/(f, t) 



1 3 2 

Vb i=l /, S =1 
2 



6 u (f,t) = — ^ qiAx^rfjiiijQiji&t) 

V " l<i<3/, 9 =l 
i<j<3 



l 3 2 

cj/(x, t) = -= V V t)z7 i g /J (f, t) 

! 2 

/o/(x, t) = -7= V (x, t), 

1 3 2 

t) = -7= Y] Y] qij(x, t)i~f 5 j iqiJ (x, t) 
1 2 

* l<i<3 f=l 
i<j<3 

where r b (6 = 1,2,3) are the Pauli matrices for the flavor indices / and g, and I and J 
denotes whether we use the local quark field (L) or the smeared field (S) to control the 
ground-state overlap. In Table II, we summarize the quantum numbers of each meson field. 
The two-point correlation functions between the interpolation fields are calculated as 

*)4(£0)> = Cu^{t), (3-31) 



x,y 
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1 3 

^( P/ (^t)p t J (y,o)) = -^c /J , 7i (t), 



x,y 



i=l 



^2(a O i(x,t)alj(y,0)) = CVj,i(t), 



x,y 



1 3 

^2{a u (x,t)a\j(y,0)) = - ^ C 7 j, 757l (t), 



x,y 



i=l 



8<j 



2<»7z(^*)»7j(j/.0)) = C /Ji75 (t) -2D /Ji75 (t), 

1 3 

J>i(*,*)4(£0)> = ^[Cun® -2D IJm (t)} 



i=l 



1 3 



x,y 



i=l 



3-32) 
3-33) 

3-34) 

3-35) 
3-36) 

3-37) 
3-38) 

3-39) 
340) 



Kj 



in terms of the connected and disconnected quark loop contributions (Tr is for the trace over 
color and spinor indices only): 



CuAt) = Yl (^i(x,t)rqi(x^{y,0)rqj(y,0) 

= ~ ^2 (Tr [Gu(x, t; y, 0)rG JI (y, 0; x , t)T]) (T = 275, ij h 1, 1757*, 17*7^), (3-41) 

£>/j,r(*) = - J] tqj&t)rqfa t) ^0)7^(y, 0) 



X ( < Tr [G n (x, t- x, t)r] - X (Tr \g h (x', t; x>, t)r 



.i.LJ 



Tr [Gjj(y, 0; y, 0)r] - ^ (Tr [g Ji7 (j?, 0; y', 0)P 

1? 

(r = ry 5 , Z7i, 1, Z7 5 7i, I7i7,-). 



(3-42) 
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Here Gu(x, t; y, t') is the propagator of the four dimensional quark field 

G<zf> a '(x,t;y,t')= [D-\x,t;y,t')]' 



Gr/' a '(x,t;y,t') = J2 FC ' C "&< 



F c c (x',y), 
D-\x',t;y,t' 
D-\x',t;y',t') 



(3-43) 
(3-44) 

(3-45) 

F c "' c 'G/,f), (3-46) 



C",C" „l „/ 



x',y 



where D 1 is written in terms of the inverse of the five dimensional matrix D D ^ F (Eq. (3-2)): 
D~\x,y) = (q{x)q{y) 

=E 



1 - 75 j. . 1 + 75 r 

0s,0 H r 8 L 1 



E 



2 

1 - 7, 



2 

1+75 



i/j(x, s)ip(y,s') 



1 + 75, . 1-75, 

s n H Os,L s -l 



A f 1 + 75 



2 

1 - 75 



8 s , Ls ^ (3-47) 



F(x,y) is the smearing function which is defined in Eq.(3-14). 

c, c', c", c w are the color indices and a, a' are the spin indices. We apply the zero- 
momentum projection to obtain the meson mass from meson energy: Eg = */ 



mi 



»2 



m meson , by summing over spatial volume x, x', y, y'. In eq. (3-41), the sum over y is stochas- 
tically evaluated by the Z 2 noise source at t — 0, while the sums over x and y in (3-42) are 
evaluated by Z 2 source spreads over the space-time volume, c.f. (3-19) and (3-20): 



N aoi 



-1 - 1 * noise 

'noise -, - - - 
-1 ^noise 

— E E^^^ )^^^ 

loise , _ _ _ 
n=l aj,3/,z 

0)^75^^ £ , i ; 0)7 5 r] ) (iV noisc — > 00) , 

J2WGu(x, t- y, 0)rG n (y, 0; x, t)r]), 
1 

-1 noise 



,r.y 



(3-48) 



■^(Tr[G 7/ (x,t;a;,t)r]) (JV n 



00 



(3-49) 
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The dagger (f ) is taken only for color and spinor (and not for space-time) indices, and we 
use the 75 hermiticity, 7 5 _D~ 1 75 = [-D" 1 ]^, of the propagator (347) in (3-48). The trace over 
color and spinor indices is exactly carried out by solving the quark propagator 3x4 times 
each for a random source. 

3.5. Meson mass fit 

Throughout this paper, we assume that the one particle state is the ground state for 
quantum numbers I and J PC , for compatibility with to the interpolation operator in Table 
II. This assumption is not entirely true for some cases. For example, a p meson may decay 
into pions. In our simulation, quarks are heavy with the lightest quark mass about half the 
strange quark mass, and confined in a relatively small (~ 2 fm) 3 box. Many of the decay 
processes would not occur in this setting since the decaying particles have energies above the 
threshold. Also we restrict ourselves to degenerate up and down quarks, Nf = 2, so that a 
meson such as ao can not decay due to exact symmetry. 

To extract the meson masses, the following two analyses are carried out. 

(A) Standard method: 

Only the ground state of mass mo is assumed to exist in the correlation function 
(OsOs), which is fitted by the hyperbolic cosine function reflecting the periodic bound- 
ary condition for a meson at t — T; 

J^(Or(x,t)O\(y,0)) = ^-\(0\ Ol \O(p = 0)>| 2 [ e -™o* + e -™oCr-*)] ? (/ = L, S) 

(3-50) 

for sufficiently large t and T — t. Although our main results will be obtained from the 
smeared-quark case, / = S, we also analyze local quark case to monitor the excited- 
state contamination. The fitting range of t is determined so that the effective meson 
mass becomes independent of the time. We also avoid a too large t for which the 
statistical error becomes large and the results become unreliable. 

(B) Variational method: 28 )' 29 ^ 

In this case, we also assume the first excited state of mass mo*. Both the local 
(/, J = L) and the smeared (/, J = S) interpolation fields are used to construct the 
correlation function (OjOj). The 2x2 matrix, 

x(t)= ( E s ,y(o L (x,t)oUy,o)) E 3 ,y(o L (x,t)oUy,o)) ) (351) 

is normalized at a reference time £ to reduce the statistical error, then is diagonalized 
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as 




X- x '\t Q )X{t)X-y 2 {t )^ l A °^ l °) u . (3 . 5 2) 



The eigenvalues are fit clS cL function of t, 



-m t , -m (T-t) . , N 

Xo (t, t ) = , +e p-— f ' <0 ^ T/2 e -o(t-to)^ (3 . 53 ) 

v ' e -m o t _|_ e -m o (T-t ) \ J y I 

Ao-(Mo) = ° t ' _„„ <^V»-"-<"> (3 . 54) 



e -m o *t -|- e -m o *(T-t ) 

to obtain the masses of the states. 

The second method, called the variational method, is employed to extract the ground- 
state energy precisely and to determine the amount of excited state contamination. 
To fit A(t, to) using eq. (3-53), without unknown amplitudes in front of the exponentials, 
to should be sufficiently large to ignore the higher excited states. By monitoring X(t, to), 
we verify, for our choice of to, that such contamination is not apparent within the 
current statistics. As an example ao case is shown in Fig. 2. to = 2 (squares) is 
chosen for the final results as X(t, to) for to = 1 (circles) can't be fit to a linear function 
of t — to meaning the meson propagator is not a single exponential, while those of 
t > 2 (diamonds, triangles) have much larger error bars. If the number of available 
configurations were larger, we would have observed the effect from the second excited 
state and should have calculated for more variations of interpolation field. This point 
may be important for future investigations with larger statistical sample. 

For another example, the effective mass of ao, which we will define in (4-1), obtained from 
the two methods is plotted in Fig. 3. The effective mass obtained from the variational method 
(squares) has the smallest statistical error, which is consistent with that of standard method 
using a smeared-smeared interpolation field (filled circles). For the variational method, the 
plateau appears after a smaller time distance when the excited-state is separated from the 
ground-state. The global fits to the plateaux are almost identical to each other. The clear 
signal of contamination from larger excited states for the local-local interpolation field (open 
circles) is observed. The identical central values and error bars from the standard single 
exponential fit and the variational method indicate that the effect from excited states is 
small for both methods with these settings. 

We analyze all masses by both methods and compare the results to estimate the sys- 
tematic uncertainty due to higher excited-states. We also explore the first excited state for 
pseudoscalar and vector mesons, 7r* and p*, using the variational method. 
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Fig. 2. i dependence of ao eigenvalue for nif = 0.02. We chose to 
the contamination from the higher excited states. 



2 (squares) by determining 



3.6. Decay constant 

The leptonic decay constant can be obtained from the amplitude of the two-point corre- 
lation function of a meson. We analyze decay constants for a pion, ir* and p mesons. Their 
respective decay constants, f n , f n * and f p can be defined through the conserved axial and 
vector currents, A b p (x) and V b (x), 

f m = (0\A b 4 (x)\O(p = 0)> = Z A (0\A b 4 (x)\O(p = 0)> (O = (3-55) 
f p m p e t = (0\V b (x)\p(p = 0)> = Z v (0\V b (x)\p(p = 0)) (i = 1,2,3) (3-56) 

where q is the polarization vector of the vector meson state, and Za and Zy are the matching 
factors between the lattice local currents, 

A\{x) =q(x)r\ l5 q(x), (3-57) 
V b (x) =q(x)T b lp q(x), (3-58) 

and an appropriate renormalization is used scheme in the continuum QCD, which, in our 
case, is MS at p, — 2 GeV. 

For /„. and f n *, the first matrix element in (3-55) can be related to pseudoscalar density 
P b (x,t) = g(x)r 6 75g(x) using the (flavor nonsinglet) AWTI, 

^(0|4(x)O(0)|0) = 2(m / + m rcs )(0|P 6 (x)O(0)|0), (3-59) 
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Fig. 3. Comparison of the effective mass of ao obtained by the two methods. Open (filled) circles 
show the results from the local-local (smeared-smeared) interpolation field in the standard 
method while squares show the effective mass obtained in the variational method. 

which leads to 

f o m 2 o = 2(m f + m ICS )(0\P b \O(p = 6)) (0 = ir,n*). (3-60) 

The actual determination of the decay constants is performed by the standard method 
(C) for a pion and p meson, and the variational method (D) for a pion and ir* meson: 

(C) Standard method 

In this case we assume the {hl^l) an d (plPl) correlation functions contain only prop- 
agation of the ground-state. (n L 7r L ) and (plPl) are fitted by a standard hyperbolic 
cosine function: 

J2Mx,t)4(y,0)) = ^-\(0\P£W= 0))| 2 [e-^ + e- m ^] 



2m„ 

8(m f + m res ) 2 

V 
2m 



[e~ m ^ + e-^V-t)] , (3-61) 



EW^AO)) = ^-\<WlW= 0))| 2 [ e -M + e -*»pCr-t)] 



x,y 



Vf^, ^_ mpt + e - mp (T-t^ ; (3 .g 2) 



to extract the quantities m„ f n , m p , and f p /Z 



v ■ 



1(3 



(D) Variational method 

In this case, the second excited state, tt*, in the correlation function of the local meson 
operator, {til^l), is also taken into account, (ttltil) is fitted by a double hyperbolic 
cosine function: 

J> L (f,t)7ri(y,0)} = _^|<0|P L >(p = 0)>| 2 [ e -^ + e -m,(T- t) ] 

+ -^H(0|P?|7r*(p = 0))| 2 \e- m ^ + e -"v( T -*)l 
2m 7r » L J 

8(m f + m ies y 1 J SK + m^) 2 L J 1 ; 

In this fitting procedure, we first determine m n and m n * by the variational method, 
(B) in the previous subsection and then fit the two-point function data to (3-63) to 
determine f n and f n * using the results from the first fitting. 

3.7. Chiral Extrapolation 

To obtain the masses and decay constants of various mesons at the physical quark mass 
point, mf = m Uj d, 25 ^ we need to extrapolate the numerical value calculated at heavier quark 
mass points. As the number of simulation points is limited and the statistical error is too 
large, we do not use the fitting formula of chiral perturbation theory at the next leading 
order or higher in this work. 

As a crude estimation of the mass of rj' at the physical point, we examine the formula 
valid in the lowest-order approximation from the flavor singlet AWTI given by eq. (3-59): 

m 2 , =C + Ci(m f + m rcs ) (AWTI type), (3-64) 

We also examine the simplest linear extrapolation for all meson masses as well as the decay 
constants, 

O = C + Ci(rrif + m TCS ) (linear type) (3-65) 
where O is either a meson mass or a decay constant. 

§4. Numerical results 

4.1. Mass of p meson and lattice scale 

First we analyze the mass of a p meson using the methods (A) and (B) and determine 
the lattice scale from m p assuming that it is a stable particle, which is true for the relatively 
heavy quark in the small box used in our simulation. In Fig. 4, the effective mas of a rho 
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Table III. m p 



rrif 


m p 


to 


^min 


^max 


method 


0.02 


0.5741(39) 




5 


n 


(A) 




0.5729(41) 


5 


fo + 1 


12 


(B) 




0.5425(64) a 




5 


16 


(A) 


0.03 


0.5979(40) 




7 


14 


(A) 




0.5984(34) 


5 


*o + l 


16 


(B) 




0.5946(58) a 




6 


16 


(A) 


0.04 


0.6385(39) 




6 


14 


(A) 




0.6379(35) 


5 


*o + l 


16 


(B) 




0.6323(70) a 




7 


16 


(A) 



( a ) These values are obtained from (I, J) = [L, W) correlation functions and quoted by. 25 ) 

meson, taken from the damping rate between meson propagators at two neighboring times, 
m c p S jj(t + 1/2), which is defined as 

(Oi(x, t)0\{y, 0)^ e -»<„(t+!) t + e -m^„(t+|) (T-t) 

is plotted in the top panels (method (A) and method (B) are shown in the left and right 
panels, respectively). The bottom panel shows an eigenvalue of the ground-state obtained 
from the variational method. The results of m p obtained from the standard hyperbolic 
cosine fit (method (A)) and the variational method (method (B)) are listed in Table III. 
The masses obtained from both methods are consistent with each other within statistical 
error for all m^, and the ground-state mass can be successfully extracted using the smeared 
operator. 

We perform linear extrapolation for both results and obtain m p at the physical quark 
mass point (rrif = m Ui( j). The result of the chiral extrapolation is shown in Fig. 5 and 
Table IV. The values obtained from both methods at the physical quark mass point are 
consistent within statistical error; we choose the value from method (B) as our main value. 
The lattice scale determined from m p =775.49 MeV 1 ) is 

a-\ = 1.537(26) GeV . (4-2) 

We have measured the potential energy between static quarks and extracted the Sommer 
scale r from the potential r /a = 4. 278(54). 25 ) Using a mp we obtained 

r P hys = 0.5491(93) fm , (4-3) 
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Fig. 4. Effective mass of p vs t using method (A) (left) and method (B) (right), and p eigenvalue 
vs t — to (bottom). Lines show the globally fitted result with errors and the ranges of t. 

which is somewhat larger than previously estimated values by ~ 10%. Although r is one of 
the most precisely determined dimensionful quantities in the lattice QCD, its experimental 
value is not known; thus, we could not judge whether our larger value is close to the physical 
value in QCD or whether it reflects some systematic errors, which we discuss in a later 
section. 

By increasing the statistical sample size, the lattice scale changed from that we reported 
in our previous paper. 25 ) Accordingly, the physical quark mass point, m/ = m Ut d, may 
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Fig. 5. m p vs rrif 
Table IV. m p at the physical quark mass point (m/ = m u d). 





< [GeV] 


a mp [fm] 


method 


0.5073(85) 


1.528(26) 


0.1291(22) 


(A) 


0.5044(85) 


1.537(26) 


0.1284(22) 


(B) 



change. However, we use the old value of m u d as the physical quark mass point in this 
paper. This is because the number of quark mass points newly obtained in this work is not 
sufficient to repeat the same analysis as before in which we used the formula of ChPT up 
to the next to the leading order. We will discuss the decay constant and the excited-state 
meson, p*, in later sections. 

4.2. Pion mass 

In Fig. 6, we plot the effective mass of a pseudoscalar meson obtained by method (A) 
on the left, that obtained by method (B) on the right, and the ground-state eigenvalue 
obtained using method (B) on the bottom panel. Table V summarizes the values of the pion 
mass obtained by both methods. By using 5-10 times more statical samples than in the 
previous analysis and extracting the ground-state information from the meson propagator 
over shorter time distance, which becomes possible using smeared operators, the statistical 
errors decrease to approximately half of those in the previous results. The fact that the 
reduction of the error size is closer to or even larger than that expected from the increase 
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Fig. 6. Pion effective mass vs. t using method (A) (left) and method (B) (right), and pion eigen- 
value vs. t — to (bottom). Lines show fit values, errors and ranges. 



in the number of statistical samples, 1/2 > 1 /\/5-10, suggests that the smearing itself does 
not necessarily cause the smaller statistical error for a pseudoscalar meson. Rather, we 
could determine the extent of the excited-state contamination using smeared operators with 
different overlaps with the states. In fact, our new results are consistent within statistical 
error with the previous results. We will discuss the decay constant and the excited-state 
meson, 7r*, later. 
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Table V. m. 





m n 


to 


^min 


^max 


method 


0.02 


0.2940(14) 




5 


14 


(A) 




0.2934(13) 


5 


*o + l 


11 


(B) 




0.2902(28) a 




9 


16 


(A) 


0.03 


0.3596(11) 




6 


16 


(A) 




0.3581(10) 


5 


to + l 


16 


(B) 




0.3575(19) a 




9 


16 


(A) 


0.04 


0.4075(11) 




7 


16 


(A) 




0.4092(11) 


5 


*o + l 


16 


(B) 




0.4094(25) a 




9 


16 


(A) 



These values are obtained from (I, J) = [L, W) correlators and quoted by. 25 ) 
4.3. Mass of do 

From experiments, there are two flavor-non-singlet scalar mesons, ao(980) and ao(1450), 
in nature. Although these are unstable particles in the more realistic Nf = 2 + 1 case, we 
assume a stable one-particle state to be the ground state in the scalar meson sector in our 
Nf = 2 case with a relatively heavy quark and small space-time, ao meson spectrum results 
previously obtained by lattice QCD calculation seem to fall roughly into two categories, 43 ^ 
studies reporting lighter masses of ~ 1 GeV 42 -*' 44 -* and those reporting heavier masses ~ 1.5 
GeV _26),43),45),46) p revious RBC results 26 ) are m ao = 1.58(34) GeV by the analysis of unitary 
points and 1.51 (19) GeV by partially quenched analysis. 

Fig. 7 (left) shows the effective mass of ao, Fig. 7 (right) shows the eigenvalue of the 
ground-state using the variational method (B). The numerical values are listed in Table VI, 
in which we also quote the previous RBC values. 26 ^ 

Our new results for the mass of ao are significantly lighter than those of previous results, 
as shown in Table VI. Since the QCD ensemble used in both investigations is the same, 
this discrepancy must originate from the difference in measuring the meson operator. In 
the previous calculation the meson interpolation field was constructed from quark fields at 
a point. Although the point operator was convenient for theoretical investigation in the 
previous study, it is not necessarily optimal for extracting the ground state. In fact, as 
shwon in the left panel of Fig. 7, the effective mass of the point operator (open symbols) is 
very large at a short distance, which implies a large amount of excited-state contamination 
in the point operator. On the other hand, the effective mass obtained using the smeared 
operator (filled symbols) reaches plateaux earlier in time and coincides to that obtained from 
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Fig. 7. Effective mass of ao vs t using method (A) (left) and method (B) (right). Lines show fitted 
values, errors, and ranges. 



Table VI. m, 



m f 




to 


^min 


^max 


method 


0.02 


0.750(28) a 




2 


6 


(A) 




0.747(28) 


2 


to + 1 


10 


(B) 




0.92(9) b 




4 


10 


exponential fit 


0.03 


0.816(17) a 




2 


6 


(A) 




0.807(14) 


2 


to + 1 


6 


(B) 




0.99(10) b 




5 


10 


exponential fit 


0.04 


0.814(19) a 




4 


9 


(A) 




0.811(17) 


4 


t + l 


10 


(B) 




0.94(5) b 




5 


12 


exponential fit 



( a ) These values are obtained by uncorrelated fitting. ^ These values are obtained from (/, J) = (L, L) 
correlators and quoted in 26). 



the ground-state eigenvalue by the variational method, shown in right panel. Note that the 
size of the statistical sample is increased by a factor of five or more in this work compared 
with that in the previous report. 

Fig. 8 shows the results of extrapolation by linear fitting and Table VII shows m ao at 
the physical quark mass point. Since both methods (A) and (B) are consistent with each 
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Fig. 8. m ao vs rrif. The asterisk on the left shows the experimental values. ^ 
Table VII. m ao at the physical quark mass point {mf = m U} d)- 





m phys [ Me y] 


m ao r 


method 


0.721(54) 


1,108(85) 


3.08(23) 


(A) 


0.723(51) 


1,111(81) 


3.10(22) 


(B) 



other we choose 

mf = 1.111(81) GeV (4-4) 

from method (B) as our final value in this work. 

To clarify the discrepancies among results obtained from lattice calculations and exper- 
iments, further investigations including the calculation multiparticle scattering states and 
the strange sea quark effects (Nf = 2 + 1) are needed. 

4.4. Mass of rj' 

Before presenting the mass spectrum results of the flavor singlet pseudoscalar meson, rj', 
we check whether the theoretical expectation discussed in §2 is realized. The ratio of the 
correlation function between disconnected quark loops, D l5 (t), to the connected correlation 
function, C 75 (i), was shown to approach unity for a large time separation, which is clearly 
different from the expectation of the linear growth in the quenched QCD case (2-9). In the 
discussion, only the pion and r( states were considered when coupling to the I(J P ) = 0(0~) 
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operator, leading to 

NrD (t) e -m v it , -m ,(r-t) , s 

In Fig. 9, the ratio extracted using the smeared operator, rj' s , is plotted. Indeed the ratio 
asymptotically approaches one for the two lighter quark masses (circles, squares), although 
it is statistically uncertain at a large time distance. However the heaviest quark mass point 
(diamonds) seems to approach to a value lower than one. The mass difference between rj' 
and 7T is smaller for the heavier quark mass, and is close to zero for the heaviest quark 
mass (as we will discuss later); thus, the ratio only approaches unity at a very large t from 
(2-6). Moreover, this deviation from the simplest theoretical explanation might be due to 
the omission of the excited states such as the ir* or 0~ + glueball state, which may play a 
more significant role in the heavier-quark-mass region. From the current results, we can not 
conclude whether the deviation from unity for the heaviest quark can be explained by the 
above-mentioned arguments or is due to other reasons, for example, insufficient sampling of 
the different topological sectors. 




Fig. 9. N f D j5 (t)/C^(t) vsi. 

We now describe the r/ spectrum obtained using methods (A) and (B). Figure 10 shows 
the effective mass (left: method (A), right: method (B)), and the ground-state eigenvalue, 
and their numerical values and their fitted ranges are given in Table VIII. We did not use 
a propagator from longer distance, where the statistics are too poor and the standard error 
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Table VIII. m rf . 



rrif 


m v > 


to 


^min 


^max 


method 


0.02 


0.477(40) 




2 


5 


(A) 




0.473(50) 


2 


*o + l 


6 


(B) 


0.03 


0.571(60) 




3 


5 


(A) 




0.600(44) 


2 


t + l 


6 


(B) 


0.04 


0.497(17) 




2 


5 


(A) 




0.492(15) 


2 


t + l 


5 


(B) 



analysis would not be reliable, although the inclusion of a few more data points does not 
change the fitted results for most of the masses. Method (B) produces flatter plateaux than 
method (A) for this meson. 

As a consistency check, we also examined the temporal exponent of the ratio (2-6) to 
extract the mass of rj'. We have evaluated the effect of the finiteness of the lattice in the 
temporal direction by using the fitting formula (2-5), and found the results to be unchanged. 
Combining the measured pion mass in Table V, the values obtained are = 0.458(58), 
0.571(48), and 0.461(15) for mf = 0.02, 0.03, and 0.04, respectively. These estimations 
are slightly smaller than the results in Table VIII. One reason for this may be that the 
time range used in fitting the ratio is too short and the pion mass is overestimated, which 
causes the estimation for the mass of rf to be smaller than its actual value. Because of this 
possibility, we won't use the results obtained from the ratio fitting in our main results. 

The mass of 7/ has only slight dependence on the quark mass, as shown in Fig. 11: all 
three masses are consistent within two to three standard deviations of statistical error. Their 
central value fluctuates nonmonotonically in quark mass order. Before being convinced of 
this nonmonotonicity, we should question the reliability of the error estimation and other 
systematic uncertainties such as insufficient sampling over the topological charge since 7/ is 
likely to depend strongly on the topological charge strongly. In our simulation we use DBW2 
gauge action to reduce the size of the residual chiral symmetry breaking, m res , sacrificing 
the configuration mobility among different topological sectors to some extent. 

Although the quark mass dependence has not been resolved sufficiently clearly, we ex- 
trapolate the measured masses by the eqs. (3-64) and (3-65) to estimate the mass of 7/ at 
the physical quark mass point. The results are shown in Fig. 11 and Table IX. The central 
values of the estimation differ from each other by 15% but are within statistical error. Our 
main estimation for the mass of 7/ at the physical quark mass point is obtained from the 
variational method (B) and chiral extrapolation using the lowest order of ChPT (3-64), and 
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Fig. 10. Effective mass of r/ vs t using method (A) (left) and method (B) (right), and rj' eigenvalue 
vs t — to (bottom). 



is given by 



m phys = 819 ( 12 7) MeV 



(4-6) 



This is the first estimation of the mass of rj' performed with the two flavors of a dynamical 
(approximately) chiral fermion, which is certainly heavier than a pion, which is thought to 
be related to the chiral U(1)a anomaly. Apart from the large statistical error and the various 
systematic errors discussed above, the main results are close to the experimentally obtained 



27 




Fig. 11. vr?, vs mj (left), and vsmj (right). The open circles and squares are extrapolated 
values at the physical quark mass point using (3-64) and (3-65), respectively. The asterisk on 
the left shows the experimental values. 1 ) 



Table IX. m„' at the physical quark mass point (rrif = m Uj( j). 



rrirf 


mP, hys [MeV] 


m v >r 


method and chiral extrapolation 


0.480(78) 


738(121) 


2.05(33) 


(A) AWTI type (3-64) 


0.487(78) 


748(120) 


2.08(33) 


(A) linear type (3-65) 


0.532(82) 


819(127) 


2.28(35) 


(B) AWTI type (3-64) 


0.560(89) 


862(130) 


2.40(36) 


(B) linear type (3-65) 



mass of rj', which suggests that further improvements can be made by especially calculation 
using an Nf = 2 + 1 ensemble. 

4.5. Mass of uj 

We also examine the flavor singlet vector meson, u, using a similar procedure to that 
for 7]'. Fig. 12 shows the effective mass of u (left: method (A), right: method (B)), which 
is also listed in Table X. We are able to extract a non-zero signal, but from a shorter time 
distance; thus there may be a significant distortion from the excited states. The results for 
the lightest point, mf = 0.02, has a particularly poor signal. 

We estimated the extrapolated mass at the physical quark mass point, as shown in Fig. 13 
and Table XI. The fitting formula used is the linear extrapolation (3-65). Since the statistical 
error for the lightest point Vfif — 0.02 is large, as mentioned above, we examine two ways of 



28 



1.2 I — 
1.1 - 

0.9 - 
| 0.8- 

<n 

■| 0.7 - 
I 0.6- 
0.5 - 
0.4 - 
0.3 - 
0.2 I — 



o 


0.02 LL 


• 


0.02 SS 


□ 


0.03 LL 


■ 


0.03 SS 


o 


0.04 LL 


♦ 


0.04 SS 



1.2 I — 
1.1 - 
1 - 
0.9 - 
1 0.8 - 

a 

<n 

■| 0.7 - 

jo 

0.5 - 

0.4 - 

0.3 - 

0.2 I — 
I) 



• m^O.02, t =2 
■ m^O.03, t =2 

♦ m 1= 0.04, t =2 




Fig. 12. Effective mass of u vs t using method (A) (left) and method (B) (right). 



Table X. m. 



m f 


rriu, 


to 


^min 


^max 


method 


0.02 


0.464(48) 




3 


6 


(A) 




0.616(124) 


2 


*o + l 


4 


(B) 


0.03 


0.636(24) 




2 


5 


(A) 




0.651(28) 


2 


to + l 


4 


(B) 


0.04 


0.717(23) 




2 


5 


(A) 




0.699(29) 


2 


to + l 


5 


(B) 



the chiral extrapolation: using all three masses or the heaviest two points. In Table XI, one 
can see that the results obtained from the "method (A) 3-masses fitting" are significantly 
different from those of the "method (A) 2-masses fitting". At the physical quark mass point, 
rriu is obtained from the "method (B) 3 masses fitting", 

m phys = 790(194) MeV . (4-7) 

Our estimation for u> is consistent with the experimental value, but with a large statistical 
error ~ 25%. 

We also calculated the propagators of the flavor singlet meson scalar, /o, using the same 
quark propagator for 7/ and u>, and found that they are too noisy to extract the spectrum 
for all values of rrif. 
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Fig. 13. vs nif. The asterisk on the left shows the experimental value. 
Table XI. Estimation of m w at the physical quark mass point {jrif — ^u^d)- 





m phys [ MeV ] 


m w r 


fit method 




0.285(80) 


439(123) 


1.22(34) 


(A) 3 masses 


fit 


0.394(117) 


605(180) 


1.68(50) 


(A) 2 masses 


fit 


0.514(126) 


790(194) 


2.20(54) 


(B) 3 masses 


fit 


0.509(141) 


782(217) 


2.18(60) 


(B) 2 masses 


fit 



4.6. Pseudovector meson (a\, b\, j\, hi) spectra 

Figs. 14, 15, 16, and 17 show the effective mass obtained using method (A) (left) and 
method (B) (right), and Tables XII, XIII-XV list the results of fits for a±, bi, fi, and 
hi, respectively. Except for the hi meson propagator at mj = 0.03, the fitting procedure 
converges. 
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Fig. 14. Effective mass of a± vs i using method (A) (left) and method (B) (right). 
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Fig. 15. Effective mass of b\ vs t using method (A) (left) and method (B) (right). 
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Fig. 16. Effective mass of j\ vs t using method (A) (left) and method (B) (right). 
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Fig. 17. Effective mass of h± vs t using method (A) (left) and method (B) (right). 
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Table XII. m 0l . 



irif 


m ai 


to 


tmin 


tmax 


method 


0.02 


0.816(12) 




4 


7 


(A) 




0.808(15) 


4 


<o + l 


9 


(B) 


0.03 


0.894(16) 




5 


8 


(A) 




0.880(11) 


4 


to + 1 


9 


(B) 


0.04 


0.898(11) 




4 


8 


(A) 




0.895(12) 


4 


to + 1 


11 


(B) 


Table XIII. m bl . 


m f 


m bl 


to 


tmin 


tmax 


method 


0.02 


0.849(17) 




4 


7 


(A) 




0.848(21) 


4 


to + 1 


8 


(B) 


0.03 


0.892(14) 




4 


8 


(A) 




0.898(15) 


4 


to + 1 


8 


(B) 


0.04 


0.889(26) 




6 


9 


(A) 




0.925(13) 


4 


to + 1 


9 


(B) 


Table XIV. m fl . 


m f 


m h 


to 


tmin 


tmax 


method 


0.02 


0.801(41) 




2 


5 


(A) 




0.798(40) 


2 


to + 1 


5 


(B) 


0.03 


0.895(34) 




2 


4 


(A) 




0.893(35) 


2 


to + 1 


4 


(B) 


0.04 


0.925(43) 




2 


4 


(A) 




0.935(41) 


2 


to + 1 


5 


(B) 



These meson masses are extrapolated linearly to the physical quark mass point, rrif — 
m Ut d, and are shown in Figs. 18-21. The numerical values are summarized in Tables XVI- 
XIX for ax, bi, fi, and hi, respectively. As the masses are independent of the method used, 
within statistical error, we choose 

mj ys = 1.140(51) GeV , (4-8) 
m Phys = L 203 (64) GeV , (4-9) 
mf?* = 1.033(137) GeV , (4-10) 
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Table XV. m hl . 



m f 




m hl 


to 


^min ^max 


method 


0.02 


0. 


818(78) 




2 4 


(A) 




0. 


814(78) 


2 


t„ + l 4 


(B) 


0.04 


0. 


834(49) 




2 4 


(A) 




0. 


832(49) 


2 


t + l 4 


(B) 



m ph ys = x 225(250) GeV (4-11) 
from method (B) as our main values. 
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Fig. 18. m ai vs rrif. The asterisk on the left shows the experimental value. ^ 
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Fig. 19. mfe 1 vs TTi j . The asterisk on the left shows the experimental value. ^ 

These numbers may be compared with the experimental results for 6 1 (1235), /ii(1170), 
ai(1260), and /i (1285), the first two of which are in good agreement with the numerical 
results. However further investigations based on realistic settings are clearly needed for 
more detailed comparisons. 
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Fig. 20. nifa vs nif. The asterisk on the left shows the experimental value. 1 ) 




Fig. 21. vs mf. The asterisk ont he left shows the experimental value. 1 ) 

4.7. Excited meson (ir* , p*) masses 

In this subsection, the second excited states of the pion and p meson are discussed. Using 
method (B), we extract the eigenvalue for the second excited state, Xo*(t),0 = it, p eq. 
(3-52), which is plotted in the right panels of Fig. 22 and 23, respectively. Although we only 
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Table XVI. m ai at the physical quark mass point (m/ = m U( i). 



m ai m* s [MeV] 


m ai r Q 


method 


0.745(26) 1,146(45) 
0.742(31) 1,140(51) 


3.19(12) 
3.17(14) 


(A) 
(B) 


Table XVII. at the physical quark mass point {rut = 


m bl mf s [MeV] 


m bl r 


method 


0.807(43) 1,241(70) 
0.783(40) 1,203(64) 


3.45(19) 
3.35(17) 


(A) 
(B) 


Table XVIII. rrif 1 at the physical quark mass point {rrif = 


m h m^ ys [MeV] 


m h r 


method 


0.689(90) 1,058(139) 
0.672(88) 1,033(137) 


2.95(39) 
2.87(38) 


(A) 
(B) 


Table XIX. m/ ll at the physical quark mass point {rrif = \ 


m hi m fti yS [MeV] 


m hl r 


method 


0.802(164) 1,233(252) 
0.797(162) 1,225(250) 


3.43(70) 
3.41(69) 


(A) 
(B) 


Table XX. 






TTlf Ul n * ^min ^max 


method 


0.02 1.215(50) 5 t - 


f 1 8 


(B) 


0.03 1.211(27) 5 t - 


f 1 8 


(B) 


0.04 1.242(26) 5 t - 


f 1 8 


(B) 



use two different operators for each meson, and Ao* (t) may have a significant contribution 
from the higher excited state, we fit \o*(t) to extract the temporal exponent, nip*, or the 
mass of the excited states using eq. (3-54). The results of the fitting are shown in Tables XX 
and XXI. We checked that the results for to = 5 and to = 6 are consistent with each other. 

We performed linear extrapolation using eq. (3-65) to the physical quark mass point, 
and found that 

m^ ys = 1.791(138) GeV , (4-12) 
mp s = 2.028(131) GeV (4-13) 
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Fig. 22. Effective mass of 7r* and eigenvalue as functions of t and t — to- 




Fig. 23. Effective mass of p* and eigenvalue as functions of t and t — to- 

(see Figs. 24 and 25, Table XXII and XXIII). These states may be interpreted as 7r(1300), 
and p(1450) or p(1700). 

4.8. Decay constants 

As the last set of numerical results, we present the leptonic decay constant in this sub- 
section. The decay constant of the ground-state pion, f n , is determined using method (C), 
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Table XXI. m p * . 



m f 


m p * 


to 


^min ^max 


method 


0.02 


1.375(43) 


5 


t + l 8 


(B) 


0.03 


1.361(25) 


5 


t + 1 8 


(B) 


0.04 


1.402(31) 


5 


t + l 8 


(B) 




Fig. 24. vn n * vs. rrif. The left most star symbols show the experimental values 1 ) in the real world. 
Table XXII. m n * at the physical quark mass point (m/ = m U} d) 

m-a* m^* ys [MeV] tii^tq method 
1.165(88) 1,791(138) 4.98(38) (B) 

Table XXIII. m p * at the physical quark mass point (nif = m Uj d) 

m p * m^ ys [MeV] m p *r method 
1.319(82) 2,028(131) 5.64(36) (B) 



and fitting the smeared two-point function to formula (3-61). We also fitted the same two- 
point functions to the double exponential formula (3-63) using the values of and m^* 
determined from the variational method (method (B)) to investigate the decay constant for 
the second excited state, J^., using method (D). 

Table XXIV shows the results for each simulated quark mass. The pion mass and decay 
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Fig. 25. m p * vs. rrif The left most star symbols show the experimental values 1 ) in the real world. 

Table XXIV. f w and fa. 



m f 


m n 


u 






^min 


^max 


method 


0.02 


0.2936(13) 


0.09561(40) 






7 


14 


(C) 




0.2934(13) (fixed) 


0.09540(43) 


1.215(50) (fixed) 


0.02244(54) 


4 


14 


(D) 




0.2938(18) a 


0.09494(62) a 






9 


16 


(C) 


0.03 


0.3598(15) 


0.10350(46) 






10 


16 


(C) 




0.3581 (10) (fixed) 


0.10370(44) 


1.211(27)(fixed) 


0.03236(65) 


4 


16 


(D) 




0.3610(18) a 


0.10253(56) a 






9 


16 


(C) 


0.04 


0.4098(12) 


0.11002(39) 






8 


16 


(C) 




0.4092(11) (fixed) 


0.10964(40) 


1.242(26) (fixed) 


0.04362(61) 


1 


16 


(D) 




0.4087(16) a 


0.11059(57) a 






9 


16 


(C) 



( a ) These values are quoted in the previous paper. 25 ) 



constants are consistent with those reported in the previous paper 25 ) within statistical error. 

Although the tt* decay constant is poorly numerically determined, an interesting theo- 
retical prediction can be made. AWTI, (3-60), for n* describes the equation for its decay 
constant, 

= 2{m > + 2 m ™ ] (0\n«') ■ (4-14) 
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Fig. 26. f w * vs rrif. 

Table XXV. f n * at the physical quark mass point (m/ = m U( j) and the chiral limit (mt = —m res ). 



m f f n * 


f^ s [MeV] 




method 


m u4 0.0013(12) 


20(19) 


0.0057(53) 


(D) 


-m rcs -0.0003(13) 


-05(20) 


-0.0013(57) 


(D) 



Table XXVI. f p /Z v . 



m f 


wi p 


f P /Zv 


^min 


^max 


method 


0.02 


0.5730(96) 


0.2011(66) 


9 


13 


(C) 


0.03 


0.6035(64) 


0.2025(50) 


10 


14 


(C) 


0.04 


0.6448(51) 


0.2164(37) 


9 


14 


(C) 



If win* is not an NG boson, so m n * remains nonzero, the right-hand side vanishes at the 
chiral limit, (wif — > — wi rcs ). This prediction was checked on a lattice QCD using Wilson 
fermions. 35 '' and f n * was consistent to be zero at the chiral limit. 

Figure 26 and Table XXV show the linear extrapolation of f n *. At the chiral limit, the 
Ti* decay constant is also consistent with the theoretical prediction, i.e. , f n * — > 0. 

Next we discuss the p meson decay constant, f p . The result of the fitting using eq. (3-62) 
is shown in Table XXVI. The mass of the p extracted by this fitting is consistent with those 
obtained from methods (A) and (B) within statistical error for all wij. 
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Then the decay constant at the physical quark mass point is obtained as 

jphys = 210(15) MeV (4-15) 

by linear extrapolation (see Table XXVII). The renormalization factor, Zy, which converts 
the lattice operator into the one in the continuum for MS at \i — 2 GeV is necessary to obtain 
a physical value for the decay constants. We use Za = 0.75734(55), which was determined 
in the previous paper, 25 ) and the relation Zy = Za, assuming the good chiral symmetry of 
the current simulation. 

0.24, , 1 , 1 , 1 , 1 , 1 , 1 




quark mass 

Fig. 27. f p /Zv vsm/. 



Table XXVII. f p at the physical quark mass point (rrif = m u ,d) 



f P /Zy 


fp 


/Phys [ MeV ] 




method 


0.1800(123) 


0.1363(94) 


210(15) 


0.0583(41) 


(C) 



4.9. Systematic uncertainties 

So far we have mainly discussed error due to the limited size of the statistical sample. Our 
numerical results were obtained only at one lattice scale, in one space-time volume, for three 
quark masses heavier than the physical values, and the strange sea quark was neglected. In 
this section, various sources of systematic errors are listed and some of their magnitudes are 
very roughly estimated to compare our results with those of experiments. 
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Approximating the continuous space-time by a discrete lattice results in a discretization 
error. Using DWF, the error starts with 0(m Tes a) + 0{a 2 A^ CT) ) . The value of m ves a 
is negligibly small in our simulation compared with the large statistical error involved 
except in the case of the pion. Our results are closer to their continuum values than 
those obtained using a Wilson-type fermion on similar lattice scale. 
For quenched DWF QCD, the physical values of /„-, fx, and fx/ fn shift by ~5%, 3%, 
and 2%, respectively, when the lattice scale changes from a" 1 = 2 GeV to continuum 
limit, 36 ) which are equal or less than current statistical error. 

Because of the limited number of quark mass points calculated in our simulation, we 
restricted ourselves to using the simplest linear chiral extrapolations (3-65) and that 
obtained from the AWTI (3-64). A more appropriate extrapolation based on a larger 
number of quark mass points is the chiral fitting formula from the (partially quenched) 
chiral perturbation theory. While the mass of rj', which is investigated as the main topic 
in this work, shows little dependence on quark mass, a more precise chiral extrapolation 
to the physical quark mass point using lighter quark masses is needed to obtain more 
reliable results. 

Although our assumption, that the ground state is a one-particle state is certainly 
wrong for some quantum numbers, some of the decay channel in nature are prohibited 
in simulations using degenerate up and down quarks with heavier mass in a relatively 
small spatial box (2 fm) 3 without a strange quark. More sophisticated investigations 
such as calculating the scattering amplitudes between multiparticles are needed to 
verify our spectrum results for the decaying meson. 

Without results obtained from a larger volume, it is difficult to estimate the finite- 
volume effect, although it might be smaller than that for baryons. 
Strange sea quark effect: The number of quark flavors that play dynamical roles in the 
rj meson may be very important as seen in the WV relation , mz, oc Nf. By increasing 
Nf from 2 to 3 by including strange quark, the WV prediction for becomes ~ 20% 
larger. A strange quark is, however, heavier than up/down quarks, and the mass of rj' 
in the Nf = 2 + 1 QCD is likely to be in between the results of Nf = 2 and 3. 
Topological charge distribution and its effects to rj' meson: In our simulation, we 
deliberately used a special gauge action, DBW2, for good chiral symmetry. However, 
the autocorrelation time of the topological charge in the simulation becomes longer. 
The samples taken in our simulation may not be sufficiently long for the reliable es- 
timation of the autocorrelation time for Qt op . The growth of the binned-jackknife 
error for (Qtop^ with increasing bin size was monitored, and we estimated the au- 
tocorrelation time of roughly ~ 300 trajectories for the rrif — 0.02 ensemble and ~ 
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Fig. 28. Circles show the measured xtop as a function of m/ 25 ^' 53 ^ while squares show values 
calculated from m„' and m n , as described in the next section. The horizontal line shows the 
value obtained from a pure SU(3) YM simulation. 4 ) The dotted line shows the prediction from 
chiral perturbation theory. 



200 trajectories for nif = 0.03,0.04. Because of the less frequent tunneling between 
different topological sectors, the charge distribution sampled in our simulation may be 
statistically skewed. In fact, ^Q top ^ = -0.7(7), 1.4(6), and 1.8(4) for m f = 0.02,0.03, 
and 0.04, respectively. Note that the central value for rrif = 0.04 is more than four 
standard deviations away from zero. It is conceivable that this poor sample of the 
topological sectors causes significant systematic errors in 7/ spectrum, particularly for 
the Tfif — 0.04 ensemble. 

Figure 28 shows the topological susceptibility, Xtop in (1*4) as a function of quark 
mass. The fact that the susceptibility for all three masses is constant within two 
standard deviations implies that the simulation points are far from the lighter-quark- 
mass region, where the susceptibility may vanish as a linear function of quark mass. 
It is also possible the tunneling between different topological sectors does not occur 
sufficiently frequently, as shown in Figure 29; thus, the estimation for the susceptibility 
has a larger systematic error. 



Of course, more reliable estimation of the magnitude of these systematic errors may be 
carried out by future simulations on a finer and larger lattice using lighter quark masses with 
the strange sea quark effect, and with a larger statistical sample size. 
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Fig. 29. History of the topological charge in the same simulation as that used for Fig. 28. 

§5. Summary and discussions 

We have measured light meson propagators in all channel (flavor nonsinglet / singlet pseu- 
doscalar, vector, scalar, pseudovector, and tensor meson it, p, a , a±, b\\ r)', u, f , fi, hi) 
and estimated the ground state meson masses and some of leptonic decay constants, as well 
as, the excited state mass, in two flavors of domain wall QCD. 

The size of the statistical sample used in the calculation is increased by five to ten 
times higher than that reported previously. 25 ) By applying the gauge-invariant Wuppertal 
smearing because of quark operators for their better overlap with the ground-state, the 
statistical error of the pion and p masses is reduced by approximately 50% and the reduction 
for rj' is more than 100%, i. e., we were only able to obtain the nonzero signal by using smeared 
field. To extract values for the meson mass and decay constant by fitting the propagators, we 
use two methods, the standard and variational methods. The results of these methods are 
consistent with each other, which indicates that excited-state contamination of the ground- 
state is controlled by the smearing. 

The systematic uncertainties discussed in the previous section are difficult to estimate; 
thus, we only quote results with statistical errors. Our results linearly extrapolated to the 
physical quark mass point are 

a'} = 1.537(26) GeV , 
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r = 0.5491(93) fm 



for quantities directly related to the lattice scale, 

f p = 210(15) MeV , 
U* =20(19) MeV 

for decay constants, and 

m ao = 1.111(81) GeV , 
rrty = 819(127) MeV , 
m u = 790(194) MeV , 
m n * = 1.791(138) GeV , 
m p * = 2.028(131) GeV , 
m ai = 1.140(51) GeV , 
m bl = 1.203(64) GeV , 
m h = 1.033(137) GeV , 
m hl = 1.225(250) GeV , 

for the mass spectrum. The lattice scale is set from m p = 775.49 MeV. 

In Fig. 30, the meson masses obtained in this work are compared with the experimental 
values. ^ Horizontal bars show the experimental values and filled circles show the simulation 
results. The error bars indicate statistical errors only. 

The decay constant of the excited pseudoscalar meson turns out to be consistent with 
zero at the chiral limit as expected: 

= -05(20) MeV . (5-1) 

The flavor singlet scalar meson, /o, was too noisy to obtain its mass in our data. 

In this paper we chose the simplest noise method, complex Z2, for evaluating the quark 
loop amplitudes. More elaborate and/or sophisticated methods 37 ) -41 ) may improve the sta- 
tistical accuracy of the calculation. 

The recalculated pion mass is consistent with the previous result, but the results for 
p and ao meson masses are significantly different. The central value of m p is 10% larger, 
i.e., a~* is 10% smaller, and the error bar is reduced by 50% compared with the previous 
results. 25 ^ Both the central value and the error bar of m ao are 25% smaller than those in 
previous results. 26 ) 
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Nf=2 DWF at a _1 =1.5 GeV on (2.1 fm) 3 
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0.2 




Fig. 30. Comparison of simulation results with experimental values 1 ) in the real world. Horizontal 
bars show the experimental values and filled circles show the simulation results. The error bars 
indicate statistical errors only. The squares show the quantities used to set the lattice spacing 
and the physical quark mass point. 



We confirm that the flavor singlet pseudoscalar meson, rj', is not an NG boson, and m v i is 
not likely to be zero at the chiral limit, which is consistent with the standard understanding 
of the axial anomaly. Assuming the WV relation (1-3) is exact at N c = 3, one can calculate 
the mass gap, ttiq, and topological susceptibility, Xto P , from our values of m v i and m n , 

Xtop(WV) = Jz-m 2 , m 2 = m*-mlpp (5-2) 

The value of Xto P (WV) for Nf = 2 are plotted in Fig. 28 as squares. The horizontal line is Xtop 
obtained from a pure SU(3) YM simulation. 4 - ) For m/ = 0.02 and 0.03, Xtop(WV) is consis- 
tent with the quenched value, while rrif — 0.04 point Xto P (WV) undershoots the line signifi- 
cantly. By linearly extrapolating to the chiral limit, we obtained = (808(129) MeV) 2 and 
Xtop(WV) = (193(15) MeV) 4 , which is consistent with the quenched value 4) (191(5) MeV) 4 . 
The agreement, which may imply only small 1/N C correction, is interesting and deserves 
further investigation in future. 

These results are susceptible to various systematic errors. First, we have only two flavors 
of dynamical quarks. The omission of the strange quark and antiquark pairs in vacuum, 
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whose mass is comparable to the dynamical scale of the QCD, may skew our results signifi- 
cantly. The limited number of quark masses, three unitary points, restricted us to examining 
only the simplest function for the quark mass dependence of the physical results. Thus, the 
chiral extrapolation has a systematic error due to the omission of curvature resulting from 
the chiral logarithms and higher order terms although many of our results show little depen- 
dence on quark mass. The ensemble was generated only on a 16 3 x 32 lattice with periodic 
boundary condition in the space directions; thus, all the meson spectrum is affected by the 
"mirror" images located ~ 2 fm away from the original image in each of the three spatial 
directions. The effects may be as large as ~ 10% for the lightest quark mass points. The 
lattice discretization error in this study is small, 0(m rcs a) + 0(a 2 AQ CD ) ~ (9(1%). The 
previous careful studies 25 -* on the scaling violation show a ~ 5% level shift for a ~ 0.1 fm lat- 
tices. The omission of the isospin violation due to the differences in quark mass and electric 
charge is likely to be negligible compared with other sources of errors, but this issue can also 
be studied nonperturbatively using lattice. 30 ^ Despite the significant statistical error and 
the various remaining systematic uncertainties, this study should serve a benchmark calcu- 
lation for the statistical features of difficult physical quantities, disconnected diagrams, and 
computational feasibility tests. The results for the mass of rj' were close to the experimental 
value, indicates that further improvements can be made particularly calculation using an 
N f = 2 + 1 DWF ensemble. 33) ' 34) 
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